
library(MASS)
library(foreign)
library(Zelig)
library(nlme)
library(sandwhich)
library(lmtest)
library(AER)
library(mgcv)
library(gee)
library(RItools)
library(xtable)
#Requires \usepackage{Dcolumn} in Latex/Sweave preamble. 
library(apsrtable)
rm(list=ls())

setwd("H:/DarkFuture/ExperimentData/DarkFutureReplicationDistribution")
ExpData <- read.dta("DarkFutureMasterR_SubjB.dta", convert.underscore=FALSE)
attach(ExpData)

#PROBIT FIGURE MAIN TEXT
#ProbReject<-zelig(OfferReject~Delta*SubjAOffer+ repetition , model="probit", data=ExpData)
ProbReject<-zelig(OfferReject~Delta*SubjAOffer+ repetition , model="probit.gee", robust=TRUE, id="ID", corstr = "exchangeable", data=ExpData)
Delta7<-setx(ProbReject,Delta=1,SubjAOffer=c(seq(0,1,by=.05)))
Delta3<-setx(ProbReject,Delta=0,SubjAOffer=c(seq(0,1,by=.05)))
s.out<-sim(ProbReject,x=Delta3, x1=Delta7)
#ProbReject<-zelig(OfferReject~Delta*SubjAOffer , model="probit", data=Rep1)
ProbReject<-zelig(OfferReject~Delta*SubjAOffer , model="probit.gee", robust=TRUE, id="ID", corstr = "exchangeable" , data=subset(ExpData, repetition==1))
Delta7<-setx(ProbReject,Delta=1,SubjAOffer=c(seq(0,1,by=.05)))
Delta3<-setx(ProbReject,Delta=0,SubjAOffer=c(seq(0,1,by=.05)))
s.out1<-sim(ProbReject,x=Delta3, x1=Delta7)
#ProbReject<-zelig(OfferReject~Delta*SubjAOffer , model="probit", data=Rep2)
ProbReject<-zelig(OfferReject~Delta*SubjAOffer , model="probit.gee", robust=TRUE, id="ID", corstr = "exchangeable" , data=subset(ExpData, repetition==2 ))
Delta7<-setx(ProbReject,Delta=1,SubjAOffer=c(seq(0,1,by=.05)))
Delta3<-setx(ProbReject,Delta=0,SubjAOffer=c(seq(0,1,by=.05)))
s.out2<-sim(ProbReject,x=Delta3, x1=Delta7)
#ProbReject<-zelig(OfferReject~Delta*SubjAOffer+ repetition , model="probit", data=ExpData_trim)
ProbReject<-zelig(OfferReject~Delta*SubjAOffer+ repetition , model="probit.gee", robust=TRUE, id="ID",corstr = "exchangeable" , data=subset(ExpData,SubjAOffer>.5 ))
Delta7<-setx(ProbReject,Delta=1,SubjAOffer=c(seq(.5,1,by=.025)))
Delta3<-setx(ProbReject,Delta=0,SubjAOffer=c(seq(.5,1,by=.025)))
s.out_trim<-sim(ProbReject,x=Delta3, x1=Delta7)


#ProbReject<-zelig(OfferReject~Delta*SubjAOffer , model="probit", data=Rep1)
ProbReject<-zelig(OfferReject~Delta*SubjAOffer , model="probit.gee", robust=TRUE, id="ID",corstr = "exchangeable" , data=subset(ExpData, repetition==1 & SubjAOffer>.5))
Delta7<-setx(ProbReject,Delta=1,SubjAOffer=c(seq(.5,1,by=.025)))
Delta3<-setx(ProbReject,Delta=0,SubjAOffer=c(seq(.5,1,by=.025)))
s.out1_trim<-sim(ProbReject,x=Delta3, x1=Delta7)
#ProbReject<-zelig(OfferReject~Delta*SubjAOffer , model="probit", data=Rep2)
ProbReject<-zelig(OfferReject~Delta*SubjAOffer , model="probit.gee", robust=TRUE, id="ID",corstr = "exchangeable" , data=subset(ExpData, repetition==2 & SubjAOffer>.5))
Delta7<-setx(ProbReject,Delta=1,SubjAOffer=c(seq(.5,1,by=.025)))
Delta3<-setx(ProbReject,Delta=0,SubjAOffer=c(seq(.5,1,by=.025)))
s.out2_trim<-sim(ProbReject,x=Delta3, x1=Delta7)

setwd("H:/DarkFuture/Paper")

#mai A numerical vector of the form c(bottom, left, top, right) which gives the margin size specified in inches. 
pdf("RejectPlot.pdf", width=8.5, height=6, onefile=FALSE, paper="special")
par(mfrow=c(2,3),mai=c(.55,.51,.42,.05), cex.lab=1.2, cex.main=1.3 , font.main=1)
plot.ci(s.out,xlab="",ylab="Probability Offer Rejected", main="Both Repetitions", col=c("dark grey","black"),ylim=c(0,1))
legend(.05, 0.2, legend = c("Delta=.7","Delta=.3"), col = c("black","dark grey"), lty = c("solid"))
plot.ci(s.out1,xlab="",ylab="", main="Repetion 1", col=c("dark grey","black"),ylim=c(0,1))
mtext("Full Sample", side=3, line=2, font=2)
plot.ci(s.out2,xlab="",ylab="", main="Repetition 2", col=c("dark grey","black"),ylim=c(0,1))

plot.ci(s.out_trim,xlab="Subject A Offer",ylab="Probability Offer Rejected", main="", col=c("dark grey","black"),ylim=c(0,1))
#legend(.2, 0.2, legend = c("Delta=.7","Delta=.3"), col = c("blue","red"), lty = c("solid"))
plot.ci(s.out1_trim,xlab="Subject A Offer",ylab="", main="", col=c("dark grey","black"),ylim=c(0,1))
mtext("Offers>.5", side=3, line=1, font=2)
plot.ci(s.out2_trim,xlab="Subject A Offer",ylab="", main="", col=c("dark grey","black"),ylim=c(0,1))
dev.off()





##############################
#Densities
##############################
All_D3<-subset(ExpData,Delta==0)
All_D7<-subset(ExpData,Delta==1)
Rep1_D3<-subset(ExpData, repetition==1 & Delta==0)
Rep1_D7<-subset(ExpData, repetition==1 & Delta==1)
Rep2_D3<-subset(ExpData, repetition==2 & Delta==0)
Rep2_D7<-subset(ExpData, repetition==2 & Delta==1)

#stderr <- function(x) sqrt(var(x)/length(x))
#abline(v=quantile(mean(Rep2_D7$SubjAOffer), c(0.025, 0.975), na.rm=TRUE), lty=2)

setwd("H:/DarkFuture/Paper")
pdf("SubjAOffers.pdf", width=10, height=5, onefile=FALSE, paper="special")
par(mfrow=c(2,2),mai=c(.7,.1,.6,.05), cex.lab=1, cex.main=1.1 , font.main=1, cex.axis=.8, yaxt="n")
plot(density(All_D3$SubjAOffer, kernel="epanechnikov"), ylim=c(0,3), main="All Repetitions", font.main=2, xlab="Offers")
lines(density(All_D7$SubjAOffer, kernel="epanechnikov"), lty=2)
abline(v=mean(All_D3$SubjAOffer))
abline(v=mean(All_D7$SubjAOffer),lty=2)
title(main="All Repetitions", font.main=2)
legend(0, 3, legend = c("Delta=.7","Delta=.3"), lty = c(2, 1 ))
mtext("Offer Distributions, Density Plots", side=3, line=2, font=2, adj=1.6)
plot(density(Rep1_D3$SubjAOffer, kernel="epanechnikov"), ylim=c(0,3), main="Repetition 1",  font.main=2, xlab="Offers", font.main=2)
lines(density(Rep1_D7$SubjAOffer, kernel="epanechnikov"), lty=2)
abline(v=mean(Rep1_D3$SubjAOffer))
abline(v=mean(Rep1_D7$SubjAOffer),lty=2)
plot(density(Rep2_D3$SubjAOffer, kernel="epanechnikov"), ylim=c(0,6), main="Repetition 2",  font.main=2, xlab="Offers")
lines(density(Rep2_D7$SubjAOffer, kernel="epanechnikov"), lty=2)
abline(v=mean(Rep2_D3$SubjAOffer))
abline(v=mean(Rep2_D7$SubjAOffer),lty=2)
dev.off()









rm(list=ls())

setwd("H:/DarkFuture/ExperimentData/DarkFutureReplicationDistribution")
ExpData <- read.dta("DarkFutureMasterR_IndividualLevel.dta", convert.underscore=FALSE,convert.factors = TRUE)
attach(ExpData)


setwd("H:/DarkFuture/Paper")
pdf("PostExpStratMethod.pdf", width=10, height=5, onefile=FALSE, paper="special")
par(mfrow=c(1,2),mai=c(.9,.1,.6,.05), cex.lab=1, cex.main=1.1 , font.main=1, cex.axis=.8, yaxt="n")
plot(density(subset(ExpData$aoffer,Delta==0), kernel="epanechnikov",from=0, to=1), ylim=c(0,4), main="Offer if in A Position", font.main=2, xlab="Offer")
lines(density(subset(ExpData$aoffer,Delta==1), kernel="epanechnikov",from=0, to=1), lty=2)
abline(v=mean(subset(ExpData$aoffer,Delta==0)))
abline(v=mean(subset(ExpData$aoffer,Delta==1)),lty=2)
#title(main="All Repetitions", font.main=2)
legend(0, 3, legend = c("Delta=.7","Delta=.3"), lty = c(2, 1 ))
mtext("Post-Experiment First Period Strategies", side=3, line=2, font=2, adj=1.9)
plot(density(subset(ExpData$BREJECT_PD1,Delta==0), kernel="epanechnikov",from=0, to=1), ylim=c(0,3), main="Minimum Offer to Accept if B", font.main=2, xlab="Minimum Required Offer")
lines(density(subset(ExpData$BREJECT_PD1,Delta==1), kernel="epanechnikov",from=0, to=1), lty=2)
abline(v=mean(subset(ExpData$BREJECT_PD1,Delta==0)))
abline(v=mean(subset(ExpData$BREJECT_PD1,Delta==1)),lty=2)
dev.off()









setwd("H:/DarkFuture/ExperimentData/DarkFutureReplicationDistribution")
ExpDataM <- read.dta("Meirowitz_DarkFutureR.dta", convert.underscore=FALSE)
attach(ExpDataM)

covar <- "SubjAOffer+ClassTreat"
ProbReject<-zelig(as.formula(paste("OfferReject ~", covar)) , model="probit.gee", robust=TRUE, id="ID", corstr = "exchangeable", data=subset(ExpDataM))
Treat<-setx(ProbReject,ClassTreat=1)
Control<-setx(ProbReject,ClassTreat=0)
Meirowitz<-sim(ProbReject,x=Control, x1=Treat, num=5000)
fd.vec<-c(mean(Meirowitz$qi$fd))
fd.up<-c(quantile(Meirowitz$qi$fd, c(0.025, 0.975), na.rm=TRUE)[2])
fd.dw<-c(quantile(Meirowitz$qi$fd, c(0.025, 0.975), na.rm=TRUE)[1])
MEffect_NoCovars<-cbind(fd.dw,fd.vec,fd.up)


MEffect_NoCovars
